# Compute distance on samples with Wolbachia
d.wolbachia <- dist(otu_wolbachia[,-1, drop=FALSE], method="euclidean")

# HCA
hca.ward <- hclust(d.wolbachia,method="ward.D2")

# 2 groups
plot(hca.ward)
rect.hclust(hca.ward,k=2)

# 3 groups
plot(hca.ward, ask=FALSE)
rect.hclust(hca.ward, k=3)

# set displaying
rh <- rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))

beg_clus <- head(cumsum(c(1, lengths(rh))), -1)
beg_clus <- beg_clus + 5
y_clus <- weighted.mean(rev(hca.ward$height)[2:3], c(4, 1))

# plot
plot(hca.ward, cex=0.7) 
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)

# save plot
setwd(path_plot)
tiff("1Ec_dendogram.tiff", units="in", width=18, height=8, res=300)
plot(hca.ward, cex=0.7, xlab="Samples", sub="") 
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)
dev.off()

png("1Ec_dendogram.png", units="in", width=18, height=8, res=300)
plot(hca.ward, cex=0.7, xlab="Samples", sub="") 
rect.hclust(hca.ward, k=3, border = c("red", "blue", "grey"))
text(x=beg_clus, y=y_clus, col=c("red","blue", "grey"), labels=c("High infection", "Low infection", "Medium infection"), font=2)
dev.off()

# split into 3 groups
hca.groups <- cutree(hca.ward, k=3)

# list of groups
print(sort(hca.groups))

# convert summary to table
summary.table <- function(df, n, var){
var <- df[df$Group==n, var]
x <- summary(var)
a <- data.frame(x=matrix(x),row.names=names(x))
colnames(a) <- n
return(a)
}

# make table with stats
wolbachia.summary <- summary.table(wolbachia.groups, 1, "Wolbachia") %>% cbind(summary.table(wolbachia.groups, 2, "Wolbachia")) %>% cbind(summary.table(wolbachia.groups, 3, "Wolbachia"))

# print
wolbachia.summary %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# print levels of Group
levels(wolbachia.groups$Group %>% as.factor)

# change levels of Group
wolbachia.groups[wolbachia.groups$Group==1, "Group"] <- "Low infection"
wolbachia.groups[wolbachia.groups$Group==2, "Group"] <- "Medium infection"
wolbachia.groups[wolbachia.groups$Group==3, "Group"] <- "High infection"

# make table with stats
wolbachia.summary.labels <- summary.table(wolbachia.groups, "Low infection", "Wolbachia") %>% cbind(summary.table(wolbachia.groups, "Medium infection", "Wolbachia")) %>% cbind(summary.table(wolbachia.groups, "High infection", "Wolbachia"))

# print
wolbachia.summary.labels %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# Compute distance on samples with Wolbachia
d.wolbachia <- dist(qpcr_wolbachia[,-1, drop=FALSE], method="euclidean")

# HCA
hca.ward <- hclust(d.wolbachia,method="ward.D2")

# 3 groups
plot(hca.ward)
rect.hclust(hca.ward, k=3)

# 2 groups
plot(hca.ward)
rect.hclust(hca.ward,k=2)

# split into 3 groups
hca.groups <- cutree(hca.ward, k=2)

# list of groups
print(sort(hca.groups))

qpcr.groups <- print(sort(hca.groups)) %>% as.data.frame()

colnames(qpcr.groups) <- "Group"
qpcr.groups$Sample <- rownames(qpcr.groups)

qpcr.groups <- qpcr.groups %>% rbind(qpcr_no_wolbachia)

# print
qpcr.groups %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

rownames(qpcr) <- qpcr$True_TubeID
qpcr <- qpcr %>% merge(qpcr.groups, by="row.names")
rownames(qpcr) <- qpcr$Sample.ID
qpcr <- qpcr[,-1]
qpcr <- qpcr %>% select(-c(Sample))

# print the qpcr table
qpcr %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# make table with stats
qpcr.summary <- summary.table(qpcr, 0, "Ratio_Target_Ref") %>% cbind(summary.table(qpcr, 1, "Ratio_Target_Ref")) %>% cbind(summary.table(qpcr, 2, "Ratio_Target_Ref"))

# print
qpcr.summary %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# print levels of Group
levels(qpcr$Group %>% as.factor)

# change levels of Group
qpcr[qpcr$Group==0, "Group"] <- "No infection"
qpcr[qpcr$Group==1, "Group"] <- "Low infection"
qpcr[qpcr$Group==2, "Group"] <- "High infection"

# make table with stats
qpcr.summary.labels <- summary.table(qpcr, "No infection", "Ratio_Target_Ref") %>% cbind(summary.table(qpcr, "Low infection", "Ratio_Target_Ref")) %>% cbind(summary.table(qpcr, "High infection", "Ratio_Target_Ref"))

# print
qpcr.summary.labels %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# print wolbachia.groups (phyloseq)
wolbachia.groups %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# print qpcr.groups (qpcr)
qpcr %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# merge 
colnames(wolbachia.groups) <- c("Sample", "Infection", "Wolbachia")
colnames(qpcr) <- c("qPCR_ID", "Sample", "qPCR_Ratio_Target_Ref", "qPCR_Infection")
all.groups <- wolbachia.groups %>% merge(qpcr, by="Sample", all.x=TRUE)

all.groups %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

# optionnal : all.groups for comparison
all.groups.comparison <- wolbachia.groups %>% merge(qpcr, by="Sample")

all.groups.comparison <- all.groups.comparison %>% select(c(Sample, qPCR_ID, Infection, Wolbachia, qPCR_Ratio_Target_Ref, qPCR_Infection))

all.groups.comparison$Infection <- all.groups.comparison$Infection %>% as.factor()

all.groups.comparison$Infection <- factor(all.groups.comparison$Infection, levels = c("Low infection", "Medium infection", "High infection"))

all.groups.comparison <- all.groups.comparison[order(all.groups.comparison$Infection),]

all.groups.comparison %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

all.groups.comparison <- all.groups.comparison %>% merge(sam %>% select(c(Sample, Strain, Species)), by="Sample")

all.groups.comparison <- all.groups.comparison[order(all.groups.comparison$Infection),]

all.groups.comparison %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

means <- aggregate(qPCR_Ratio_Target_Ref ~ Infection , all.groups.comparison, median)
means$qPCR_Ratio_Target_Ref <- means$qPCR_Ratio_Target_Ref %>% round(2)

ggqqplot(all.groups.comparison$Wolbachia)

shapiro.test(all.groups.comparison$Wolbachia)

my_comparisons <- list( c("Low infection", "Medium infection"), 
                        c("Medium infection", "High infection"), 
                        c("Low infection", "High infection") )

p1 <- ggplot(all.groups.comparison, aes(x=Infection, y=qPCR_Ratio_Target_Ref, fill=Infection)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
        geom_jitter(position=position_jitter(0.5), aes(color=Infection), alpha=0.6) +
        theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")), 
       y="Ratio Target Ref (qPCR)", 
       fill=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")))+
  guides(colour=FALSE)+
  #stat_summary(fun.y = "median", geom = "point", shape = 23, size = 3, fill = "white")+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.signif", 
                     vjust=1.5,
                     label.y = c(10, 11, 12))+
    stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.format", 
                     label.y = c(10, 11, 12))+
  geom_text(data = means, aes(label = qPCR_Ratio_Target_Ref, y = qPCR_Ratio_Target_Ref + 0.6))

p1

p9 <- ggplot(all.groups.comparison, aes(x=Infection, y=Wolbachia, fill=Infection)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
  geom_jitter(position=position_jitter(0.5), aes(color=Infection), alpha=0.6) +
  theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x=expression(paste("\n", italic("Wolbachia"), "infection groups (from HCA)")), 
       y=expression(paste("\n", italic("Wolbachia"), "infection (relative abundance)")), 
       fill="Infection")+
  guides(colour=FALSE)+ 
  scale_y_continuous(breaks = seq(0, 1, by = 0.1))

p9

setwd(path_tsv)
write.table(all.groups.comparison, "1Ec_qpcr_comparison.tsv", sep="\t", row.names = FALSE, col.names=TRUE)

setwd(path_plot)
tiff("1Ec_boxplot.tiff", units="in", width=10, height=5, res=300)
p1
dev.off()

png("1Ec_boxplot.png", units="in", width=10, height=5, res=300)
p1
dev.off()

# load phyloseq object after decontam
setwd(path_rdata)
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")

# extract metadata
sam <- data.frame(ps.filter@sam_data)

# add groups in sam table
sam <- sam %>% merge(all.groups %>% select(c(Sample, Infection, Wolbachia, qPCR_Infection, qPCR_Ratio_Target_Ref)), by="Sample", all.x = TRUE)
rownames(sam)<- sam$Sample

# change levels for plot
# levels(sam$Strain) <- c(levels(sam$Strain), "Laboratory - Lavar", "Laboratory - Slab TC (Wolbachia -)")
# sam$Strain[sam$Strain=="Laboratory - Lavar"] <- "Laboratory - Lavar"
# sam$Strain[sam$Strain=="Wolbachia -"] <- "Laboratory - Slab TC (Wolbachia -)"
# sam$Strain <- droplevels(sam$Strain)
# levels(sam$Strain)

# add a column for qPCR measure
sam$qPCR_measure <- "Yes"
sam[is.na(sam$qPCR_Ratio_Target_Ref), "qPCR_measure"] <- "No"
sam[sam$qPCR_measure=="Yes" & sam$qPCR_Ratio_Target_Ref==0, "qPCR_measure"] <- "= 0 (Wolbachia -)"
sam[sam$qPCR_measure=="Yes" & sam$qPCR_Ratio_Target_Ref>0, "qPCR_measure"] <- "> 0 (Wolbachia +)"

# update sam table of phyloseq object
sample_data(ps.filter) <- sam

# transformation percent
ps_percent <- transform_sample_counts(ps.filter, function(x) x/sum(x)*100)

# select whole
ps_percent_whole <- subset_samples(ps_percent, Organ == "Whole")

# all info
p1 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Infection", shape=NA) +
  labs(x="NMDS1", y = "NMDS2")+
  stat_ellipse(geom = "polygon", level=0.95, alpha = 0.25, aes(fill = Species, color=NA))+
  scale_fill_manual(values=c("orange", "yellow","green"), labels = c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
  scale_color_manual(values=c("blue", "grey","red"), breaks=c("Low infection", "Medium infection", "High infection"))+
  scale_shape_manual(values=c(19,17,15,4,7), labels=c("Field - Bosc", "Field - Camping Europe", "Field - Guadeloupe", "Laboratory - Lavar", expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)"))))+
  theme_gray()+
  labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  guides(fill = guide_legend(order = 1),
         colour = guide_legend(order = 2),
         shape = guide_legend(order = 3))+
  labs(fill="Species", col="Wolbachia Group")+
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  theme(legend.text.align = 0)

p1$layers[[1]] <- NULL
p1[["guides"]][["shape"]][["title"]] <- "Strain"
p1[["guides"]][["colour"]][["title"]] <- expression(paste(italic("Wolbachia"), "Infection"))

p1+
  geom_point(aes(shape = Strain), size = 2)

# Wolbachia infection

p2 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Infection", shape="qPCR_measure") +
  
  scale_color_manual(values=c("blue", "#E4A7F5","red"), 
                     breaks=c("Low infection", "Medium infection", "High infection"))+
  scale_fill_manual(expression(paste(italic("Wolbachia"), "Infection")), 
                    values=c("blue", "#E4A7F5","red"), 
                    breaks=c("Low infection", "Medium infection", "High infection"))+
  scale_shape_manual("qPCR", 
                     values=c(7,2,19), 
                     labels=c(expression(paste("= 0 (", italic("Wolbachia"), "-)")), 
                              expression(paste("> 0 (", italic("Wolbachia"), "+)")), 
                              "No qPCR"))+
  labs(x="NMDS1", 
       y = "NMDS2",
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  
  stat_ellipse(geom = "polygon", 
               alpha = 0.25, 
               aes(fill=Infection, group=Infection))+
  
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         shape = guide_legend(order = 2),
         colour=FALSE)+
  
  theme(legend.text.align = 0)

p2

# species ("orange", "yellow", "green")
p3 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Species") +
  
  geom_point(size = 2) +
  
  scale_fill_manual(values=c("#F57462", "#498CF5","#559E36"), 
                    labels = c(expression(italic("Aedes aegypti")), 
                               expression(italic("Culex pipiens")), 
                               expression(italic("Culex quinquefasciatus"))))+
  
  scale_color_manual(values=c("#F57462", "#498CF5","#559E36"))+
  
  labs(x="NMDS1", 
       y = "NMDS2",
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  
  stat_ellipse(geom = "polygon", 
               alpha = 0.25, 
               aes(fill = Species))+
  
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         colour = FALSE,
         shape = guide_legend(order = 3))+
  
  
  theme(legend.text.align = 0)

p3

# Strain
new_location <- c("Field - Guadeloupe", "Laboratory - Slab TC (Wolbachia -)", 
                  "Field - Bosc", "Field - Camping Europe", "Laboratory - Lavar")
ps_percent_whole@sam_data$Strain <- factor(ps_percent_whole@sam_data$Strain, levels = new_location)

p4 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Strain") +
  
  geom_point(size = 2) +
  
  scale_fill_manual("Strain", 
                    values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
                    labels = c("Field - Guadeloupe", 
                               expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
                               "Field - Bosc", 
                               "Field - Camping Europe", 
                               "Laboratory - Lavar"))+
  scale_color_manual("Strain", 
                     values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
                     labels = c("Field - Guadeloupe", 
                                expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
                                "Field - Bosc", 
                                "Field - Camping Europe", 
                                "Laboratory - Lavar"))+
  
  labs(x="NMDS1", 
       y = "NMDS2",
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  stat_ellipse(geom = "polygon", alpha = 0.25, aes(fill=Strain, group=Strain))+
  
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         colour = FALSE,
         shape = guide_legend(order = 3))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  theme(legend.text.align = 0)

p4

# field
p5 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Field") +
  geom_point(size = 2) +
  scale_fill_manual(values=c("brown", "purple" ))+
  scale_color_manual(values=c("brown", "purple"))+
  
  labs(x="NMDS1", 
       y = "NMDS2", 
       caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  
  stat_ellipse(geom = "polygon", alpha = 0.25, aes(fill=Field, group=Field))+
  
  #coord_cartesian(xlim = c(-2,2.8), ylim=c(-2,1.5))+
  coord_cartesian(xlim = c(-2.8,2.8), ylim=c(-2.5,2.2))+
  #expand_limits(x=c(-2.5, 1.5), y=c(-2, 2.5))+
  
  guides(fill = guide_legend(order = 1),
         colour = FALSE,
         shape = guide_legend(order = 3))+
  
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)

p5

# species + strains
p6 <- plot_ordination(ps_percent_whole, bray.culex_whole, color="Strain", shape=NA) +
  labs(x="NMDS1", y = "NMDS2")+
  scale_color_manual("Strain", 
                     values=c("#00BF7D", "#5B6BF7", "#A3A500", "#F8766D", "#00B0F6"),
                     labels = c("Field - Guadeloupe", 
                                expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), "-)")),
                                "Field - Bosc", 
                                "Field - Camping Europe", 
                                "Laboratory - Lavar"))+
  scale_shape_manual(values = c(15,16,17), labels=c(expression(italic("Aedes aegypti")), expression(italic("Culex pipiens")), expression(italic("Culex quinquefasciatus"))))+
  theme_gray()+
  labs(caption=paste("stress:", bray.culex_whole$stress %>% round(2), "\n", "ellipse level: 0.95"))+
  guides(fill = guide_legend(order = 1),
         colour = guide_legend(order = 2),
         shape = guide_legend(order = 3))+
  geom_text(mapping = aes(label = Sample), size = 1.5, vjust = 3)+
  theme(legend.text.align = 0)+
  geom_point(aes(shape = Species_italic), size = 2)

p6$layers[[1]] <- NULL
p6[["guides"]][["shape"]][["title"]] <- "Species"

p6

sam2 <- sam[sam$Species=="Culex pipiens",]

means <- aggregate(Wolbachia ~ Strain , sam2, median)
means$Wolbachia <- means$Wolbachia %>% round(2)

ggqqplot(sam2$Wolbachia)

shapiro.test(sam2$Wolbachia)

my_comparisons <- list( c("Field - Bosc", "Field - Camping Europe"), 
                        c("Field - Camping Europe", "Laboratory - Lavar"), 
                        c("Field - Bosc", "Laboratory - Lavar") )

p7 <- ggplot(sam2, aes(x=Strain, y=Wolbachia, fill=Strain)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
  geom_jitter(position=position_jitter(0.5), aes(color=Strain), alpha=0.6) +
  theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x="\nStrain", 
       y=expression(paste("\n", italic("Wolbachia"), "infection (relative abundance)")),
       fill="Strain")+
  guides(colour=FALSE)+
  #stat_summary(fun = "mean", geom = "point", shape = 23, size = 3, fill = "white")+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.signif", 
                     vjust=1.5,
                     label.y = c(1.1, 1.2, 1.3))+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.format", 
                     label.y = c(1.1, 1.2, 1.3))+
  geom_text(data = means, aes(label = Wolbachia, y = Wolbachia - 0.05))

p7

means <- aggregate(qPCR_Ratio_Target_Ref ~ Strain , sam2, median)
means$qPCR_Ratio_Target_Ref<- means$qPCR_Ratio_Target_Ref %>% round(2)

wilcox.test(qPCR_Ratio_Target_Ref ~ Strain, sam2[sam2$Strain!="Field - Camping Europe",])

sam3 <- na.omit(sam2)
ggqqplot(sam3$qPCR_Ratio_Target_Ref)

shapiro.test(sam3$qPCR_Ratio_Target_Ref)

p8 <- ggplot(na.omit(sam2), aes(x=Strain, y=qPCR_Ratio_Target_Ref, fill=Strain)) + 
  geom_boxplot(alpha=0.6, outlier.shape = NA)+
  geom_jitter(position=position_jitter(0.5), aes(color=Strain), alpha=0.6) +
  theme_bw(base_size = 14)+
  scale_fill_manual(values = c("blue", "grey", "red"))+
  scale_color_manual(values = c("blue", "grey", "red"))+
  labs(x="\nStrain", y="Ratio Target Ref (qPCR)", fill="Strain")+
  guides(colour=FALSE)+
  #stat_summary(fun = "median", geom = "point", shape = 23, size = 3, fill = "white")+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.signif", 
                     vjust=1.5,
                     label.y = c(10, 11, 12))+
  stat_compare_means(comparisons = my_comparisons,
                     method="wilcox.test", 
                     label="p.format", 
                     label.y = c(10, 11, 12))+
  geom_text(data = means, aes(label = qPCR_Ratio_Target_Ref, y = qPCR_Ratio_Target_Ref + 0.4))

p8